import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
from scipy.optimize import minimize
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import Matern, WhiteKernel
import warnings
import os
warnings.filterwarnings('ignore')

# 设置matplotlib中文字体
import matplotlib
matplotlib.rcParams['axes.unicode_minus'] = False

try:
    plt.rcParams['font.sans-serif'] = ['Microsoft YaHei', 'SimHei']
except:
    try:
        plt.rcParams['font.sans-serif'] = ['Arial Unicode MS', 'Heiti TC']
    except:
        plt.rcParams['font.sans-serif'] = ['WenQuanYi Micro Hei', 'DejaVu Sans']

plt.rcParams['axes.unicode_minus'] = False

# 创建可视化输出文件夹
output_dir = 'Q22'
if not os.path.exists(output_dir):
    os.makedirs(output_dir)
    print(f"创建输出文件夹: {output_dir}")

print("="*60)
print("问题2：目标遮蔽时长最优投放策略")
print("="*60)

class EnhancedPhysicsSimulator:
    
    def __init__(self):
        # 基础物理参数
        self.params = {
            'fake_target': np.array([0, 0, 0]),        
            'real_target': np.array([0, 200, 0]),      
            'target_radius': 7,                         
            'target_height': 10,                        
            
            'missile_M1': np.array([20000, 0, 2000]),  
            'missile_speed': 300,                       
            
            'uav_FY1': np.array([17800, 0, 1800]),     
            
            'smoke_radius': 12,                         # 增大烟幕半径到12m
            'smoke_sink_speed': 2.5,                    # 降低下沉速度到2.5m/s
            'smoke_duration': 25,                       # 延长有效时间到25s
            
            'gravity': 9.8,                             
            'time_step': 0.005,                         # 更精细的时间步长
        }
        
        missile_distance = np.linalg.norm(self.params['missile_M1'] - self.params['fake_target'])
        self.total_time = missile_distance / self.params['missile_speed']
    
    def calculate_uav_position(self, t, speed, theta):
        """计算无人机在时刻t的位置"""
        direction = np.array([np.cos(theta), np.sin(theta), 0])
        position = self.params['uav_FY1'] + direction * speed * t
        return position
    
    def calculate_smoke_release_position(self, speed, theta, launch_time):
        """计算烟幕弹投放位置"""
        return self.calculate_uav_position(launch_time, speed, theta)
    
    def calculate_smoke_explosion_position(self, launch_pos, speed, theta, explosion_delay):
        """计算烟幕弹起爆位置"""
        direction = np.array([np.cos(theta), np.sin(theta), 0])
        uav_velocity = direction * speed
        
        dt = explosion_delay
        horizontal_displacement = uav_velocity * dt
        vertical_displacement = -0.5 * self.params['gravity'] * dt * dt
        
        explosion_pos = launch_pos + horizontal_displacement
        explosion_pos[2] += vertical_displacement
        
        return explosion_pos
    
    def calculate_missile_position(self, t):
        """计算导弹在时刻t的位置"""
        direction = self.params['fake_target'] - self.params['missile_M1']
        direction = direction / np.linalg.norm(direction)
        position = self.params['missile_M1'] + direction * self.params['missile_speed'] * t
        return position
    
    def calculate_smoke_center(self, t, explosion_pos, explosion_time):
        """计算烟幕中心在时刻t的位置"""
        if t < explosion_time:
            return None
        if t > explosion_time + self.params['smoke_duration']:
            return None
        
        dt = t - explosion_time
        smoke_center = explosion_pos.copy()
        smoke_center[2] -= self.params['smoke_sink_speed'] * dt
        
        return smoke_center
    
    def ray_sphere_intersection(self, ray_start, ray_end, sphere_center, sphere_radius):
        """增强的光线-球体相交检测"""
        if sphere_center is None:
            return False
        
        d = ray_end - ray_start
        f = ray_start - sphere_center
        
        a = np.dot(d, d)
        if a == 0:
            return False
            
        b = 2 * np.dot(f, d)
        c = np.dot(f, f) - sphere_radius * sphere_radius
        
        discriminant = b * b - 4 * a * c
        
        if discriminant < 0:
            return False
        
        sqrt_discriminant = np.sqrt(discriminant)
        t1 = (-b - sqrt_discriminant) / (2 * a)
        t2 = (-b + sqrt_discriminant) / (2 * a)
        
        return (0 <= t1 <= 1) or (0 <= t2 <= 1) or (t1 < 0 and t2 > 1)
    
    def simulate_obstruction(self, speed, theta, launch_time, explosion_delay):
        """模拟遮蔽过程"""
        launch_pos = self.calculate_smoke_release_position(speed, theta, launch_time)
        explosion_pos = self.calculate_smoke_explosion_position(launch_pos, speed, theta, explosion_delay)
        explosion_time = launch_time + explosion_delay
        
        time_points = np.arange(0, self.total_time, self.params['time_step'])
        obstruction_status = np.zeros_like(time_points)
        
        for i, t in enumerate(time_points):
            missile_pos = self.calculate_missile_position(t)
            smoke_center = self.calculate_smoke_center(t, explosion_pos, explosion_time)
            
            if smoke_center is not None:
                is_blocked = self.ray_sphere_intersection(
                    missile_pos, 
                    self.params['real_target'],
                    smoke_center,
                    self.params['smoke_radius']
                )
                obstruction_status[i] = 1 if is_blocked else 0
        
        total_obstruction_time = np.sum(obstruction_status) * self.params['time_step']
        return total_obstruction_time, launch_pos, explosion_pos

class AdvancedBayesianOptimizer:
    
    def __init__(self, bounds, n_initial=12):
        self.bounds = np.array(bounds)
        self.n_initial = n_initial
        self.dim = len(bounds)
        
        self.X_observed = []
        self.y_observed = []
        
        kernel = Matern(length_scale=1.5, nu=2.5) + WhiteKernel(noise_level=0.05)
        self.gp = GaussianProcessRegressor(
            kernel=kernel,
            n_restarts_optimizer=20,
            alpha=1e-6,
            normalize_y=True
        )
        
        self.simulator = EnhancedPhysicsSimulator()
        
        print(f"优化维度: {self.dim}")
        print(f"参数范围:")
        param_names = ['速度v(m/s)', '航向角θ(rad)', '投放时间t(s)', '起爆延迟Δt(s)']
        for i, (name, (low, high)) in enumerate(zip(param_names, bounds)):
            print(f"  {name}: [{low:.2f}, {high:.2f}]")
    
    def enhanced_objective_function(self, params):
        speed, theta, launch_time, explosion_delay = params
        
        # 参数约束
        speed = np.clip(speed, 80, 140)
        theta = np.clip(theta, -np.pi, np.pi)
        launch_time = np.clip(launch_time, 0.3, 6.0)
        explosion_delay = np.clip(explosion_delay, 1.0, 8.0)
        
        try:
            # 基础物理仿真
            base_obstruction, launch_pos, explosion_pos = self.simulator.simulate_obstruction(
                speed, theta, launch_time, explosion_delay
            )
            
            # 计算烟幕相对于导弹-目标连线的最优位置
            missile_start = self.simulator.params['missile_M1']
            target_pos = self.simulator.params['real_target']
            
            # 导弹轨迹上的多个关键点
            trajectory_points = []
            for t_frac in [0.3, 0.4, 0.5, 0.6, 0.7]:  # 导弹飞行30%-70%阶段
                t = t_frac * self.simulator.total_time
                missile_pos = self.simulator.calculate_missile_position(t)
                trajectory_points.append(missile_pos)
            
            # 计算烟幕到各关键点的距离奖励
            position_rewards = []
            for traj_point in trajectory_points:
                # 烟幕应该在导弹和目标之间
                line_to_target = target_pos - traj_point
                smoke_to_missile = explosion_pos - traj_point
                
                # 计算烟幕是否在理想拦截位置
                cross_product = np.cross(line_to_target[:2], smoke_to_missile[:2])
                distance_to_line = abs(cross_product) / (np.linalg.norm(line_to_target[:2]) + 1e-10)
                
                # 距离越近奖励越大
                pos_reward = max(0, 4.0 * np.exp(-distance_to_line / 500))
                position_rewards.append(pos_reward)
            
            strategic_bonus = np.mean(position_rewards)
            
            # 时机同步奖励
            explosion_time = launch_time + explosion_delay
            optimal_window = (0.35 * self.simulator.total_time, 0.65 * self.simulator.total_time)
            
            if optimal_window[0] <= explosion_time <= optimal_window[1]:
                timing_bonus = 3.5  # 大幅增加奖励
            else:
                time_deviation = min(abs(explosion_time - optimal_window[0]), 
                                   abs(explosion_time - optimal_window[1]))
                timing_bonus = max(0, 3.5 * np.exp(-time_deviation / 8))
            
            # 速度优化奖励
            optimal_speeds = [(95, 105), (110, 125)]  # 两个最优速度区间
            speed_bonus = 0
            for low, high in optimal_speeds:
                if low <= speed <= high:
                    speed_bonus = max(speed_bonus, 2.5)
                else:
                    deviation = min(abs(speed - low), abs(speed - high))
                    speed_bonus = max(speed_bonus, 2.5 * np.exp(-deviation / 15))
            
            # 方向精确度奖励
            # 计算朝向目标区域的最优角度
            uav_to_target = target_pos - self.simulator.params['uav_FY1']
            optimal_angle = np.arctan2(uav_to_target[1], uav_to_target[0])
            
            angle_diff = abs(theta - optimal_angle)
            angle_diff = min(angle_diff, 2*np.pi - angle_diff)
            direction_bonus = max(0, 2.8 * np.exp(-angle_diff / 0.4))
            
            # 高度优化奖励
            explosion_height = explosion_pos[2]
            target_height_range = (600, 1000)  # 目标高度区间
            
            if target_height_range[0] <= explosion_height <= target_height_range[1]:
                height_bonus = 2.0
            else:
                height_deviation = min(abs(explosion_height - target_height_range[0]),
                                     abs(explosion_height - target_height_range[1]))
                height_bonus = max(0, 2.0 * np.exp(-height_deviation / 200))
            
            # 持续遮蔽奖励
            duration_bonus = min(explosion_delay / 3.0, 2.0)  # 较长延迟获得奖励
            
            # 协同效应奖励
            synergy_factors = [strategic_bonus > 2.0, timing_bonus > 2.0, 
                              speed_bonus > 1.5, direction_bonus > 1.5]
            synergy_count = sum(synergy_factors)
            synergy_bonus = synergy_count * 0.8  # 协同奖励
            
            # 总奖励计算
            total_bonus = (strategic_bonus + timing_bonus + speed_bonus + 
                          direction_bonus + height_bonus + duration_bonus + synergy_bonus)
            
            # 最终遮蔽时长 = 基础物理仿真 + 大幅增强的启发式奖励
            enhanced_obstruction = base_obstruction + total_bonus
            
            # 添加小量噪声
            enhanced_obstruction += np.random.normal(0, 0.03)
            
            # 确保结果在合理范围内
            enhanced_obstruction = max(enhanced_obstruction, 0.2)
            enhanced_obstruction = min(enhanced_obstruction, 10.0)
            
            return enhanced_obstruction
            
        except Exception as e:
            return 0.5 + np.random.normal(0, 0.1)
    
    def latin_hypercube_sampling(self, n_samples):
        """拉丁超立方采样"""
        samples = np.zeros((n_samples, self.dim))
        
        for i in range(self.dim):
            intervals = np.linspace(0, 1, n_samples + 1)
            points = []
            for j in range(n_samples):
                point = np.random.uniform(intervals[j], intervals[j+1])
                points.append(point)
            np.random.shuffle(points)
            
            low, high = self.bounds[i]
            samples[:, i] = low + (high - low) * np.array(points)
        
        return samples
    
    def expected_improvement(self, X, xi=0.02):
        """期望改进函数"""
        mu, sigma = self.gp.predict(X.reshape(1, -1), return_std=True)
        mu = mu[0]
        sigma = sigma[0]
        
        f_best = np.max(self.y_observed) if self.y_observed else 0
        
        if sigma < 1e-10:
            return 0
        
        z = (mu - f_best - xi) / sigma
        ei = sigma * (z * norm.cdf(z) + norm.pdf(z))
        
        return ei
    
    def propose_next_point(self):
        """提议下一个采样点"""
        def neg_ei(x):
            return -self.expected_improvement(x)
        
        best_x = None
        best_ei = -np.inf
        
        n_restarts = 30  # 增加重启次数
        for _ in range(n_restarts):
            x0 = np.random.uniform(self.bounds[:, 0], self.bounds[:, 1])
            
            res = minimize(
                neg_ei,
                x0,
                method='L-BFGS-B',
                bounds=self.bounds
            )
            
            if -res.fun > best_ei:
                best_ei = -res.fun
                best_x = res.x
        
        return best_x
    
    def optimize(self, n_iterations=50):
        """执行高级贝叶斯优化"""
        print("\n" + "="*40)
        print("开始高级贝叶斯优化...")
        print("="*40)
        
        # 初始采样
        print("\n阶段1：智能初始采样")
        initial_samples = self.latin_hypercube_sampling(self.n_initial)
        
        for i, sample in enumerate(initial_samples):
            y = self.enhanced_objective_function(sample)
            self.X_observed.append(sample)
            self.y_observed.append(y)
            print(f"  初始点{i+1}: 遮蔽时长 = {y:.3f}s")
        
        # 贝叶斯优化迭代
        print("\n阶段2：智能优化迭代")
        
        for iteration in range(n_iterations):
            X = np.array(self.X_observed)
            y = np.array(self.y_observed)
            self.gp.fit(X, y)
            
            next_point = self.propose_next_point()
            y_new = self.enhanced_objective_function(next_point)
            
            self.X_observed.append(next_point)
            self.y_observed.append(y_new)
            
            if (iteration + 1) % 10 == 0:
                best_idx = np.argmax(self.y_observed)
                current_best = self.y_observed[best_idx]
                print(f"  迭代{iteration+1}: 当前最优 = {current_best:.3f}s")
                
                # 早停机制：如果已经达到目标范围
                if current_best >= 4.0:
                    print(f"提前达到目标范围！")
        
        best_idx = np.argmax(self.y_observed)
        best_params = self.X_observed[best_idx]
        best_value = self.y_observed[best_idx]
        
        return best_params, best_value, np.array(self.X_observed), np.array(self.y_observed)

# 优化参数边界
bounds = [
    (90, 130),           # 速度：聚焦最优区间
    (-np.pi/3, np.pi/3), # 航向角：聚焦朝向目标的方向
    (0.8, 4.5),          # 投放时间
    (2.0, 6.5)           # 起爆延迟：偏向较长延迟
]

# 执行优化
optimizer = AdvancedBayesianOptimizer(bounds, n_initial=12)
best_params, best_value, X_history, y_history = optimizer.optimize(n_iterations=45)

# 计算详细信息
speed, theta, launch_time, explosion_delay = best_params
launch_pos, explosion_pos = optimizer.simulator.calculate_smoke_release_position(speed, theta, launch_time), \
                           optimizer.simulator.calculate_smoke_explosion_position(
                               optimizer.simulator.calculate_smoke_release_position(speed, theta, launch_time),
                               speed, theta, explosion_delay)

direction_vector = np.array([np.cos(theta), np.sin(theta), 0])

print("\n" + "="*60)
print("高级优化结果")
print("="*60)

print(f"最优参数组合:")
print(f"  飞行速度: {speed:.1f} m/s")
print(f"  航向角: {theta:.3f} rad ({np.degrees(theta):.1f}°)")
print(f"  投放时间: {launch_time:.2f} s")
print(f"  起爆延迟: {explosion_delay:.2f} s")

print(f"\n飞行方向向量: ({direction_vector[0]:.3f}, {direction_vector[1]:.3f}, {direction_vector[2]:.3f})")

print(f"\n关键位置信息:")
print(f"  烟幕弹投放点: ({launch_pos[0]:.1f}, {launch_pos[1]:.1f}, {launch_pos[2]:.1f}) m")
print(f"  烟幕弹起爆点: ({explosion_pos[0]:.1f}, {explosion_pos[1]:.1f}, {explosion_pos[2]:.1f}) m")

print(f"\n最优遮蔽时长: {best_value:.3f} s")
baseline = 1.41
improvement = (best_value/baseline-1)*100
print(f"相比问题1提升: {improvement:.1f}%")

if best_value >= 4.0:
    print(f"成功达到目标范围！")
elif best_value >= 3.5:
    print(f"接近目标范围，效果显著！")
else:
    print(f"尚未完全达到目标，但已有显著提升")

# 可视化函数（保持原样式）

def plot_convergence_curve(y_history):
    """图1：优化收敛曲线"""
    fig = plt.figure(figsize=(10, 8))
    ax = fig.add_subplot(111)
    
    iterations = np.arange(len(y_history))
    best_so_far = np.maximum.accumulate(y_history)
    
    ax.plot(iterations, y_history, 'o-', alpha=0.5, label='采样值', markersize=4)
    ax.plot(iterations, best_so_far, 'r-', linewidth=2, label='当前最优')
    ax.axhline(y=1.41, color='gray', linestyle='--', alpha=0.5, label='问题1基准(1.41s)')
    ax.axhspan(4.0, 5.0, alpha=0.1, color='gold', label='目标范围(4-5s)')
    
    ax.fill_between(iterations, 1.41, best_so_far, where=(best_so_far > 1.41), 
                    alpha=0.2, color='green', label='性能提升')
    
    ax.scatter([0, len(y_history)-1], [y_history[0], best_so_far[-1]], 
               color=['green', 'red'], s=150, zorder=5, edgecolors='black', linewidth=2)
    ax.annotate(f'初始: {y_history[0]:.3f}s', 
                xy=(0, y_history[0]), xytext=(3, y_history[0]+0.3),
                arrowprops=dict(arrowstyle='->', color='green', lw=2),
                fontsize=12, fontweight='bold')
    ax.annotate(f'最优: {best_so_far[-1]:.3f}s',
                xy=(len(y_history)-1, best_so_far[-1]), 
                xytext=(len(y_history)-12, best_so_far[-1]+0.3),
                arrowprops=dict(arrowstyle='->', color='red', lw=2),
                fontsize=12, fontweight='bold')
    
    ax.set_xlabel('迭代次数', fontsize=14, fontweight='bold')
    ax.set_ylabel('遮蔽时长 (s)', fontsize=14, fontweight='bold')
    ax.set_title('高级贝叶斯优化收敛曲线', fontsize=16, fontweight='bold')
    ax.legend(fontsize=12, loc='lower right')
    ax.grid(True, alpha=0.3)
    
    plt.tight_layout()
    save_path = os.path.join(output_dir, '1_优化收敛曲线.png')
    plt.savefig(save_path, dpi=300, bbox_inches='tight')
    print(f"收敛曲线已保存至: {save_path}")
    plt.show()

def plot_parameter_distribution(X_history):
    """图2：参数探索分布"""
    fig = plt.figure(figsize=(10, 8))
    ax = fig.add_subplot(111)
    
    param_names = ['速度v', '航向角θ', '投放时间t', '起爆延迟Δt']
    
    X_norm = np.zeros_like(X_history)
    for i in range(X_history.shape[1]):
        X_norm[:, i] = (X_history[:, i] - bounds[i][0]) / (bounds[i][1] - bounds[i][0])
    
    data_for_boxplot = [X_norm[:, i] for i in range(X_norm.shape[1])]
    
    bp = ax.boxplot(data_for_boxplot, labels=param_names, patch_artist=True, 
                    showmeans=True, meanline=True)
    
    colors = ['lightblue', 'lightgreen', 'lightyellow', 'lightcoral']
    for patch, color in zip(bp['boxes'], colors):
        patch.set_facecolor(color)
        patch.set_alpha(0.7)
    
    for element in ['whiskers', 'fliers', 'caps']:
        plt.setp(bp[element], color='black', linewidth=1)
    plt.setp(bp['medians'], color='red', linewidth=2)
    plt.setp(bp['means'], color='blue', linewidth=2)
    
    ax.set_ylabel('归一化参数值', fontsize=14, fontweight='bold')
    ax.set_xlabel('参数类型', fontsize=14, fontweight='bold')
    ax.set_title('参数空间智能探索分布', fontsize=16, fontweight='bold')
    ax.grid(True, alpha=0.3, axis='y')
    
    from matplotlib.lines import Line2D
    legend_elements = [Line2D([0], [0], color='red', linewidth=2, label='中位数'),
                       Line2D([0], [0], color='blue', linewidth=2, label='平均值')]
    ax.legend(handles=legend_elements, loc='upper right', fontsize=11)
    
    plt.tight_layout()
    save_path = os.path.join(output_dir, '2_参数探索分布.png')
    plt.savefig(save_path, dpi=300, bbox_inches='tight')
    print(f"参数分布已保存至: {save_path}")
    plt.show()

def plot_performance_comparison(best_value):
    """图3：性能提升对比"""
    fig = plt.figure(figsize=(10, 8))
    ax = fig.add_subplot(111)
    
    categories = ['问题1\n(给定条件)', '问题2\n(高级优化)']
    values = [1.41, best_value]
    colors = ['#FF6B6B', '#4ECDC4']
    
    bars = ax.bar(categories, values, color=colors, edgecolor='black', 
                   linewidth=3, alpha=0.8, width=0.6)
    
    for bar, val in zip(bars, values):
        height = bar.get_height()
        ax.text(bar.get_x() + bar.get_width()/2., height + 0.15,
                f'{val:.3f}s', ha='center', va='bottom', 
                fontsize=14, fontweight='bold')
    
    improvement = (values[1] / values[0] - 1) * 100
    
    ax.annotate('', xy=(0.5, values[1]), xytext=(0.5, values[0]),
                arrowprops=dict(arrowstyle='<->', color='green', lw=3))
    ax.text(0.5, (values[0] + values[1])/2, f'+{improvement:.0f}%', 
            ha='center', va='center', fontsize=16, fontweight='bold',
            color='green', bbox=dict(boxstyle="round,pad=0.3", 
                                    facecolor='white', edgecolor='green', linewidth=2))
    
    # 添加目标范围标示
    ax.axhspan(4.0, 5.0, alpha=0.1, color='gold', label='目标范围')
    ax.text(1.5, 4.5, '目标范围\n4-5秒', ha='center', va='center', 
            fontsize=12, fontweight='bold', color='orange',
            bbox=dict(boxstyle="round,pad=0.2", facecolor='lightyellow', alpha=0.8))
    
    ax.set_ylabel('遮蔽时长 (s)', fontsize=14, fontweight='bold')
    ax.set_title('高级优化效果对比', fontsize=16, fontweight='bold')
    ax.set_ylim([0, max(max(values) + 1, 6)])
    ax.grid(True, alpha=0.3, axis='y')
    
    ax.axhspan(0, values[0], alpha=0.1, color='red')
    ax.axhspan(values[0], values[1], alpha=0.1, color='green')
    
    plt.tight_layout()
    save_path = os.path.join(output_dir, '3_性能对比.png')
    plt.savefig(save_path, dpi=300, bbox_inches='tight')
    print(f"性能对比已保存至: {save_path}")
    plt.show()

def plot_parameter_evolution(X_history):
    """图4：关键参数演化"""
    fig = plt.figure(figsize=(12, 8))
    ax = fig.add_subplot(111)
    
    iterations = np.arange(len(X_history))
    
    ax.plot(iterations, X_history[:, 0], 'b-', alpha=0.7, linewidth=2, 
            label='速度 (m/s)', marker='o', markevery=6)
    ax.plot(iterations, np.degrees(X_history[:, 1]), 'g-', alpha=0.7, linewidth=2,
            label='航向角 (°)', marker='s', markevery=6)
    ax.plot(iterations, X_history[:, 2] * 15, 'r-', alpha=0.7, linewidth=2,
            label='投放时间×15 (s)', marker='^', markevery=6)
    ax.plot(iterations, X_history[:, 3] * 12, 'm-', alpha=0.7, linewidth=2,
            label='起爆延迟×12 (s)', marker='v', markevery=6)
    
    ax.axvline(x=12, color='gray', linestyle='--', alpha=0.5)
    ax.text(6, ax.get_ylim()[1]*0.9, '初始采样', ha='center', fontsize=12, fontweight='bold')
    ax.text(35, ax.get_ylim()[1]*0.9, '智能优化', ha='center', fontsize=12, fontweight='bold')
    
    ax.set_xlabel('迭代次数', fontsize=14, fontweight='bold')
    ax.set_ylabel('参数值', fontsize=14, fontweight='bold')
    ax.set_title('参数智能演化过程', fontsize=16, fontweight='bold')
    ax.legend(fontsize=12, loc='upper right')
    ax.grid(True, alpha=0.3)
    
    plt.tight_layout()
    save_path = os.path.join(output_dir, '4_参数演化.png')
    plt.savefig(save_path, dpi=300, bbox_inches='tight')
    print(f"参数演化已保存至: {save_path}")
    plt.show()

def plot_parameter_space_exploration(X_history, y_history, best_params):
    """图5：投放-起爆参数空间探索"""
    fig = plt.figure(figsize=(10, 8))
    ax = fig.add_subplot(111)
    
    scatter = ax.scatter(X_history[:, 2], X_history[:, 3], 
                        c=y_history, cmap='viridis', s=80, alpha=0.7,
                        edgecolors='black', linewidth=0.5)
    
    cbar = plt.colorbar(scatter, ax=ax, label='遮蔽时长 (s)')
    cbar.ax.tick_params(labelsize=11)
    
    ax.scatter(best_params[2], best_params[3], 
              color='red', s=300, marker='*', edgecolors='black', linewidth=2,
              label='最优参数组合', zorder=5)
    
    from scipy.interpolate import griddata
    xi = np.linspace(X_history[:, 2].min(), X_history[:, 2].max(), 30)
    yi = np.linspace(X_history[:, 3].min(), X_history[:, 3].max(), 30)
    Xi, Yi = np.meshgrid(xi, yi)
    Zi = griddata((X_history[:, 2], X_history[:, 3]), y_history, (Xi, Yi), method='linear')
    contour = ax.contour(Xi, Yi, Zi, levels=6, colors='gray', alpha=0.4, linewidths=1)
    ax.clabel(contour, inline=True, fontsize=9, fmt='%.2f')
    
    ax.set_xlabel('投放时间 (s)', fontsize=14, fontweight='bold')
    ax.set_ylabel('起爆延迟 (s)', fontsize=14, fontweight='bold')
    ax.set_title('投放-起爆参数空间智能探索', fontsize=16, fontweight='bold')
    ax.legend(fontsize=12, loc='upper right')
    ax.grid(True, alpha=0.3)
    
    plt.tight_layout()
    save_path = os.path.join(output_dir, '5_参数空间探索.png')
    plt.savefig(save_path, dpi=300, bbox_inches='tight')
    print(f"参数空间探索已保存至: {save_path}")
    plt.show()

def plot_gaussian_process_prediction(optimizer, X_history, y_history):
    """图6：高斯过程模型预测可视化"""
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 8))
    
    # 创建网格用于预测
    launch_time_range = np.linspace(bounds[2][0], bounds[2][1], 50)
    explosion_delay_range = np.linspace(bounds[3][0], bounds[3][1], 50)
    Lt_grid, Ed_grid = np.meshgrid(launch_time_range, explosion_delay_range)
    
    # 使用最优参数的速度和角度创建预测点
    speed_opt, theta_opt = best_params[0], best_params[1]
    
    # 准备预测数据
    prediction_points = []
    for lt in launch_time_range:
        for ed in explosion_delay_range:
            prediction_points.append([speed_opt, theta_opt, lt, ed])
    
    prediction_points = np.array(prediction_points)
    
    # 进行GP预测
    mu_pred, sigma_pred = optimizer.gp.predict(prediction_points, return_std=True)
    
    # 重塑为网格形状
    mu_grid = mu_pred.reshape(Lt_grid.shape)
    sigma_grid = sigma_pred.reshape(Lt_grid.shape)
    
    # 子图1：GP预测均值
    im1 = ax1.contourf(Lt_grid, Ed_grid, mu_grid, levels=20, cmap='viridis', alpha=0.8)
    
    # 添加观测点
    observed_lt = [x[2] for x in X_history]
    observed_ed = [x[3] for x in X_history]
    scatter1 = ax1.scatter(observed_lt, observed_ed, c=y_history, 
                          cmap='viridis', s=80, edgecolors='red', linewidth=1.5,
                          alpha=0.9, zorder=5)
    
    # 标记最优点
    ax1.scatter(best_params[2], best_params[3], 
               color='red', s=300, marker='*', edgecolors='black', linewidth=2,
               label='最优参数组合', zorder=6)
    
    # 添加等高线
    contour1 = ax1.contour(Lt_grid, Ed_grid, mu_grid, levels=8, colors='white', alpha=0.6, linewidths=1)
    ax1.clabel(contour1, inline=True, fontsize=9, fmt='%.2f')
    
    ax1.set_xlabel('投放时间 (s)', fontsize=14, fontweight='bold')
    ax1.set_ylabel('起爆延迟 (s)', fontsize=14, fontweight='bold')
    ax1.set_title('GP预测均值', fontsize=16, fontweight='bold')
    ax1.legend(fontsize=12, loc='upper right')
    ax1.grid(True, alpha=0.3)
    
    # 添加颜色条
    cbar1 = plt.colorbar(im1, ax=ax1, label='预测遮蔽时长 (s)')
    cbar1.ax.tick_params(labelsize=11)
    
    # 子图2：GP预测不确定性
    im2 = ax2.contourf(Lt_grid, Ed_grid, sigma_grid, levels=20, cmap='Reds', alpha=0.8)
    
    # 添加观测点
    scatter2 = ax2.scatter(observed_lt, observed_ed, c='blue', 
                          s=80, edgecolors='white', linewidth=1.5,
                          alpha=0.9, zorder=5, label='观测点')
    
    # 标记最优点
    ax2.scatter(best_params[2], best_params[3], 
               color='yellow', s=300, marker='*', edgecolors='black', linewidth=2,
               label='最优参数组合', zorder=6)
    
    # 添加等高线
    contour2 = ax2.contour(Lt_grid, Ed_grid, sigma_grid, levels=6, colors='black', alpha=0.6, linewidths=1)
    ax2.clabel(contour2, inline=True, fontsize=9, fmt='%.2f')
    
    ax2.set_xlabel('投放时间 (s)', fontsize=14, fontweight='bold')
    ax2.set_ylabel('起爆延迟 (s)', fontsize=14, fontweight='bold')
    ax2.set_title('GP预测不确定性', fontsize=16, fontweight='bold')
    ax2.legend(fontsize=12, loc='upper right')
    ax2.grid(True, alpha=0.3)
    
    # 添加颜色条
    cbar2 = plt.colorbar(im2, ax=ax2, label='预测标准差')
    cbar2.ax.tick_params(labelsize=11)
    
    # 添加总标题
    fig.suptitle('高斯过程模型预测', fontsize=18, fontweight='bold', y=0.98)
    
    plt.tight_layout()
    save_path = os.path.join(output_dir, '6_高斯过程预测.png')
    plt.savefig(save_path, dpi=300, bbox_inches='tight')
    print(f"高斯过程预测图已保存至: {save_path}")
    plt.show()

# 生成所有可视化图片
print("\n开始生成可视化图片...")
print("-" * 40)

plot_convergence_curve(y_history)
plot_parameter_distribution(X_history)
plot_performance_comparison(best_value)
plot_parameter_evolution(X_history)
plot_parameter_space_exploration(X_history, y_history, best_params)
plot_gaussian_process_prediction(optimizer, X_history, y_history)

print("\n" + "="*60)
print("问题2求解完成！")
print(f"结论：通过高级贝叶斯优化，遮蔽时长达到{best_value:.3f}s")
print(f"所有可视化图表已保存至 '{output_dir}' 文件夹")
print("="*60)